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ULTRASOUND THERAPY 



FIELD OF THE INVENTION 

The invention relates to medical systems and, more particularly, to non-invasive 
5 application of focused ultrasound energy to subjects such as humans, and in particular to 
the brain of a human subject. 

BACKGROUND OF THE INVENTION 

Treatment of tissues lying at specific locations within the skull may be limited to 

10 removal or ablation. While these treatments have proven effective for certain localized 
disorders, such as tumors, they involve delicate, time-consuming procedures that may 
result in destruction of otherwise healthy tissues. These treatments are generally not 
appropriate for disorders in which diseased tissue is integrated into healthy tissue, except 
in instances where destruction of the healthy tissue will not unduly effect neurologic 

15 function. 

The noninvasive nature of ultrasound surgery has special appeal in the brain 
where it is often desirable to destroy or treat deep tissue volumes without disturbing 

healthy tissues. Focused ultrasound beams have been used for noninvasive surgery in 

f 

many other parts of the body. Ultrasound penetrates well through soft tissues and, due to 
20 the short wavelengths (1.5 mm at 1 MHz), it can be focused to spots with dimensions of a 
few millimeters. By heating, e.g., using ultrasound, tumorous or cancerous tissue in the 
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abdomen, for example, it is possible to ablate the diseased portions without significant 
damage to surrounding healthy tissue. 

SUMMARY OF THE INVENTION 
5 In general, in one aspect, the invention provides a method of delivering ultrasound 

signals. The method includes providing an image of at least a portion of a subject 
intended to receive ultrasound signals between sources of the ultrasound signals and a 
desired region of the subject for receiving focused ultrasound signals, identifying, from 
the image, physical characteristics of different layers of material between the sources and 

10 the desired region, and determining at least one of phase corrections and amplitude 

corrections for the sources depending on respective thicknesses of portions of each of the 
layers disposed between each source and the desired region. 

Implementations of the invention may include one or more of the following 
features. The physical characteristics are associated with material type and at least one of 

1 5 material density and material structure, the identifying further comprising identifying 
thicknesses of the layers. The phase corrections are determined in accordance with 
propagation characteristics of each of the layers. The propagation characteristics are 
determined based upon the material type and at least one of the material density and the 
material structure of each of the respective layers. The layers are identified using values 
20 associated with portions of the image. The values are intensities of the portions of the 
image. The phase corrections are determined using a three-layer model of a skull of the 
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subject Two of the three layers are assumed to have approximately identical speeds of 
sound, Cj, therein, with the other layer having a speed of sound c\\ therein, wherein the 
phase corrections are determined using a phase shift determined according to: 

5 where Cn is a speed of sound in the n^ layer, and D„ is a thickness of the n* layer, and 
wherein the speeds of sound in the layers are determined according to: 

c 0 c u lovjJJ 

where di, cfc, d 3 , are thicknesses of the three layers, $(p) is a measured phase shift as a 
function of density, and p is density. 

10 Further, implementations of the invention may include one or more of the 

following features. The physical characteristics are associated with x-ray attenuation 
coefficients, \x„ The material between the sources and the desired region is bone. The 
phase corrections are related to the attenuation coefficient by a phase function including 
parameters derived at least partially experimentally. Each phase correction equals M + 

15 BE(l/ji(x)) + CE(l/^x)) 2 , where ^i(x) is the attenuation coefficient as a function of 
distance x along a line of propagation between each source and the desired region, and 
where M, B, and C are derived at least partially experimentally. The amplitude 
corrections are related to the attenuation coefficient by an amplitude function including 
parameters derived at least partially experimentally. Each amplitude correction is related 

20 to N + FEn(x) + G£fa(x)) 2 , where |i(x) is the attenuation coefficient as a function of 

3 
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distance x along a line of propagation between each source and the desired region, and 

where N, F, and G are derived at least partially experimentally. 

Further, implementations of the invention may include one or more of the 

following features. The layers are identified according to both material density and 
5 material structure. Providing the image includes producing the image using magnetic 

resonance imaging. Providing the image includes producing the image using computer 

tomography. The sources are piezoelectric transducer elements. Both phase and 

amplitude corrections are determined. 

In general, in another aspect, the invention provides a system for delivering 
10 ultrasound signals. The system includes an apparatus configured to analyze an image of 

at least a portion of a subject intended to receive ultrasound signals between sources of 

the ultrasound signals and a desired region of the subject for receiving focused 

ultrasound signals, the apparatus configured to determine, from the image, information 

about different layers of the at least a portion of the subject, and an array of sources of 
15 ultrasound signals having at least one of their relative phases and their amplitudes set in 

accordance with the information about each layer of the at least a portion of the subject 

provided by the apparatus. 

Implementations of the invention may include one or more of the following 

features. The phases are set in accordance with propagation characteristics of each layer 
20 of the at least a portion of the subject. The propagation characteristics are dependent 

upon the material type and at least one of the material density and the material structure 
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of each layer of the at least a portion of the subject. The apparatus is configured to 
identify the layers using values associated with portions of the image. The values are 
intensities of the portions of the image. The apparatus is configured to determine the 
information about different layers of bone. The apparatus is configured to determine the 
5 phase corrections using a three-layer model of a skull of the subject. The information is 
associated with an x-ray attenuation coefficient, \i: The phase corrections are related to 
the attenuation coefficient by a phase function including parameters derived at least 
partially experimentally. The amplitude corrections are related to the attenuation 
coefficient by an amplitude function including parameters derived at least partially 

10 experimentally. 

Further, implementations of the invention may include one or more of the 
following features. The system further includes a magnetic resonance imager coupled to 
the apparatus and configured to produce the image. The system further includes a 
computer tomography imager coupled to the apparatus and configured to produce the 

15 image. The sources are piezoelectric transducer elements. 

In general, in another aspect, the invention provides a computer program product 
residing on a computer readable medium and comprising instructions for causing a 
computer to analyze an image of at least a portion of a subject to receive ultrasound 
signals between sources of the ultrasound signals and a desired region of the subject for 

20 receiving focused ultrasound signals to identify, from the image, physical characteristics 
of layers of material between the sources and the desired region, and to determine at least 
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one of phase corrections and amplitude corrections for the sources depending on 
respective thicknesses of portions of each of the layers disposed between each source and 
the desired region. 

Implementations of the invention may include one or more of the following 
5 features. The phase corrections are determined in accordance with propagation 

characteristics of each of the layers. The propagation characteristics are dependent upon 
the material type and at least one of the material density and the material structure of each 
of the respective layers. The layers are identified according to both material density and 
material structure. The computer program product further includes instructions for 

10 causing a computer to produce the image using magnetic resonance imaging. The 
computer program product further includes instructions for causing a computer to 
produce the image using computer tomography. The instructions for causing a computer 
to identify layers of materials are for causing the computer to identify the layers of 
materials based upon intensities of portions of the image. 

15 Further, implementations of the invention may include one or more of the 

following features. The layers are identified using values associated with portions of the 
image. The values are intensities of the portions of the image. The layers analyzed are 
layers of bone. The phase corrections are determined using a three-layer model of a skull 
of the subject. Two of the three layers are assumed to have approximately the same speed 

20 of sound, q, therein, with the other layer having a speed of sound c H therein, wherein the 
phase corrections are determined using a phase shift determined according to: 
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where c n is a speed of sound in the layer, and D n is a thickness of the layer, and 
wherein the speeds of sound in the layers are determined according to: 

c 0 c„ 360JD 

5 where dj , d 2 , d 3 , are thicknesses of the three layers, ^(p) is a measured phase shift as a 
function of density, and p is density. 

Further, implementations of the invention may include one or more of the 
following features. The physical characteristics are associated with x-ray attenuation 
coefficients, n. The phase corrections are related to the attenuation coefficient by a phase 
10 function including parameters derived at least partially experimentally. Each phase 

correction equals M + BS(l/ji(x)) + CS(l/fi(x)) 2 , where p(x) is the attenuation coefficient 
as a function of distance x along a line of propagation between each source and the 
desired region, and where M, B, and C are derived at least partially experimentally. The 
amplitude corrections are related to the attenuation coefficient by an amplitude function 
15 including parameters derived at least partially experimentally. Each amplitude correction 
is related toN + FE]Lt(x) 4- GS(|i(x)) 2 , where n(x) is the attenuation coefficient as a 
function of distance x along a line of propagation between each source and the desired 
region, and where N, F, and G are derived at least partially experimentally. 

In general, in another aspect, the invention provides a method of providing 
20 ultrasound signals into a subject from at least one source of an array of sources of 
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ultrasound signals. The method includes (a) transmitting ultrasound energy of a selected 
frequency from a selected source into the subject, (b) receiving superimposed reflections 
of the transmitted energy, the reflections being from an outer surface of the subject and at 
least one interface inside the subject, (c) repeating (a) and (b) using ultrasound energy of 
5 frequencies other than the selected frequency, (d) determining a frequency difference 
between frequencies associated with relative extrema of the received reflections, and (e) 
using the determined frequency difference and a thickness, of at least a portion of 
material between the selected source and a desired region in the subject for receiving 
focused ultrasound energy signals, to determine a phase correction for the selected source. 

10 Implementations of the invention may include one or more of the following 

features. The method further includes (f) providing an image of at least a portion of a 
subject intended to receive ultrasound energy signals between sources of the energy 
signals and the desired region, and (g) identifying, from the image, the thickness of at 
least a portion of material between the selected source and the desired region. The 

15 method further includes repeating (a) - (e) for each of the sources other than the selected 
source. The phase correction is determined according to: 

A* = 27rf[(d/co).(l/(2Af))] 
where A$ is the phase correction, f is a frequency to be transmitted, d is the thickness, c 0 
is the speed of sound in water, and Af is the frequency difference between like extrema. 

20 In general, in another aspect, the invention provides logic for use in a system for 

providing ultrasound energy into a living subject from an array of sources of ultrasound 
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energy signals. The logic is configured to control apparatus to (a) transmit ultrasound 
energy of a selected frequency from a selected source into the subject, (b) receive 
superimposed reflections of the transmitted energy, the reflections being from an outer 
surface of the subject and at least one interface inside the subject, (c) repeat (a) and (b) 
5 using ultrasound energy of frequencies other than the selected frequency, (d) determine a 
frequency difference between frequencies associated with relative extrema of the received 
reflections, and (e) use the determined frequency difference and a thickness, of at least a 
portion of material between the selected source and a desired region in the subject for 
receiving focused ultrasound energy signals, to determine a phase correction for the 

10 selected source. 

Implementations of the invention may include one or more of the following 
features. The logic is further configured to cause the apparatus to (f) provide an image of 
at least a portion of a subject intended to receive ultrasound energy signals between 
sources of the energy signals and the desired region, and (g) identify, from the image, the 

15 thickness of at least a portion of material between the selected source and the desired 
region. The logic is further configured to cause the apparatus to repeat (a) - (e) for each 
of the sources other than the selected source. The logic is configured to cause the 
apparatus to determine the phase correction according to: 
A<|)==27tf [(d/c 0 )-(l/(2Af))] 

20 where A<|> is the phase correction, f is a frequency to be transmitted, d is the thickness, cq 
is the speed of sound in water, and Af is the frequency difference between like extrema. 
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In general, in another aspect, the invention provides a method of delivering 
ultrasound signals. The method includes providing an image of at least a portion of a 
subject intended to receive ultrasound signals between sources of the ultrasound signals 
and a desired region of the subject for receiving focused ultrasound signals, identifying, 
5 from the image, physical characteristics of portions of material between the sources and 
the desired region, and using a spectral propagation analysis to determine at least one of 
phase corrections and amplitude corrections for the sources depending on the respective 
physical characteristics of the portions of the material disposed between each source and 
the desired region. 

10 Implementations of the invention may include one or more of the following 

features. The physical characteristics include density and thickness. The identifying 
identifies the physical characteristics of different layers of the portions of the material. 
The physical characteristics include angle of a layer interface relative to a reference. 

Implementations of the invention may also include one or more of the following 

15 features. The spectral propagation analysis includes a wavevector-frequency domain 
projection algorith^n. The using includes comparing projected phases from the sources to 
the desired region accounting for the material and ignoring the material. The using 
assumes perpendicular incidence of ultrasound. The method further includes determining 
an incident angle of ultrasound to the material. The identifying identifies the physical 

20 characteristics of different layers of the portions of the material and the determining an 
incident angle determines incident angles for each layer of the material. 
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In general, in another aspect, the invention provides a system for delivering 
ultrasound signals. The system includes an apparatus configured to analyze an image of 
at least a portion of a subject intended to receive ultrasound signals between sources of 
the ultrasound signals and a desired region of the subject for receiving focused ultrasound 
5 signals, the apparatus configured to determine, from the image using a spectral projection 
analysis, information about the at least a portion of the subject, and an array of sources of 
ultrasound signals having at least respective ones of their relative phases and their 
amplitudes set in accordance with the information about the at least a portion of the 
subject provided by the apparatus. 

10 Implementations of the invention may include one or more of the following 

features. The apparatus comprises means for identifying physical characteristics of 
multiple layers of the portion of the subject, the physical characteristics including density 
and thickness. The physical characteristics further include incident angles of the layers 
relative to respective sources, the apparatus further comprising means for comparing 

15 projected phases of the signals from the sources to the desired region accounting for the 
portion of the subject and ignoring the portion of the subject. The apparatus comprises 
means for performing a wavevector-frequency domain projection algorithm. 

In general, in another aspect, the invention provides a computer program product 
residing on a computer readable medium and including computer-readable, computer 

20 executable instructions for causing a computer to analyze an image of at least a portion of 
a subject to receive ultrasound signals between sources of the ultrasound signals and a 
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desired region of the subject for receiving focused ultrasound signals to identify, from the 
image, physical characteristics of respective portions of material between the sources and 
the desired region, and to determine, using a wavevector-frequency domain projection 
algorithm, at least one of phase corrections and amplitude corrections for the sources 
5 depending on respective physical characteristics of portions of the material disposed 
between each source and the desired region. 

Implementations of the invention may include one or more of the following 
features. The physical characteristics include density and thickness. The instructions for 
causing the computer to analyze cause the computer to identify the physical 

10 characteristics of different layers of the portions of the material. The physical 
characteristics include angles of layer interfaces relative to the sources. 

Implementations of the invention may also include one or more of the following 
features. The instructions for causing the computer to analyze cause the computer to 
compare projected phases from the sources to the desired region accounting for the 

15 portions of the material and ignoring the portions of the material. The computer program 
product further includes instructions for causing the computer to determine the projected 
phases accounting for the portions of the material assuming perpendicular incidence of 
ultrasound. The computer program product further includes instructions for causing the 
computer to determine the projected phases accounting for the portions of the material by 

20 determining and using incident angles of ultrasound relative to the layers of the portions 
of the material. The instructions for causing the computer to determine incident angles 
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cause the computer to determine incident angles for each layer of the portions of the 
material. 

Various aspects of the invention may provide one or more of the following 
advantages. Ultrasound can be focused accurately within an intact skull, e.g., for 
5 ultrasound therapy. Different skulls, e.g., different skull thicknesses, densities, and/or 
structures, can be accommodated for ultrasound therapy. Real-time adjustments to 
ultrasound therapy can be made. Effects on phase and/or amplitude of energy passing 
through bone (or other tissue) may be determined and used to compensate the phase 
and/or amplitude of energy applied to the bone (or other tissue). Incident angles of 

10 ultrasound relative to a skull can be accounted for and excitation characteristics of 

ultrasound transducers adjusted accordingly. A wavevector-frequency domain projection 
algorithm can be used to determine ultrasound transducer excitation characteristics for 
providing ultrasound non-invasively into a skull. Correction factors for ultrasound array 
transducers can be determined in an efficient, timely manner. 

15 These and other advantages of the invention, along with the invention itself, will 

be more fully understood after a review of the following figures, detailed description, and 
claims. 

BRIEF DESCRIPTION OF THE FIGURES 
20 FIG. 1 is a schematic diagram of an ultrasound therapy system according to the 

invention. 
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FIG. 2 is a 3-dimensional rendering of a portion of a patient's skull. 
FIG. 3 is a flow diagram of a process of determining excitation correction factors 
and exciting transducer elements using the determined factors. 

FIG. 4 is a schematic diagram of another ultrasound therapy system according to 
5 the invention. 

FIG. 5 is a flow diagram of a process of obtaining phase shift factors using the 
system shown in FIG. 4. 

FIG. 6 is schematic diagram showing planar field rotation in a wavevector- 
frequency domain, defining a new projection plane. 
10 FIG. 7 is a schematic ray diagram of ultrasound propagation through non-parallel 

layers of different-density materials. 
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DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 

Energy, such as ultrasound energy, may be propagated into a subject such as a 
human or animal, for various therapeutic reasons. Arrays of radiating transducer 
elements may be used to transmit energy into the subject, and the amplitudes of the 
5 signals transmitted by the elements can affect how much energy penetrates the subject 
and the relative phases and amplitudes of energies transmitted can help to focus a 
distribution of energy in the subject. It has been found that the phases and amplitudes of 
ultrasound signals may be affected/distorted by many causes such as different material 
properties of materials through which the signals are propagated. For example, different 

10 types of materials within a patient, e.g., bone, muscle, and fat, have different propagation 
and attenuation constants of ultrasound energy. Also, other material properties, such as 
the density and structure of materials in a patient, may affect the phase and/or attenuation 
of ultrasound signals. Because the signals from different transducer elements may 
encounter different thicknesses and contours of materials and possibly air-filled or liquid- 

15 filled pockets between transducer elements and a region to be imaged/treated, the phases 
of the signals from the transducer elements will often be distorted. The resulting phases 
and/or amplitudes will therefore often be different than desired if the transmission phases 
and/or amplitudes are not compensated for this distortion. 

Embodiments of the invention provide techniques for compensating for phase 

20 distortions and attenuation variances when treating a patient using ultrasound. Using 
imaging techniques such as computer tomography (CT), magnetic resonance imaging 
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(MRI), etc., properties such as thickness, density, and structure of materials are 
determined. It has been discovered that each of these properties affects the phase 
distortion of materials, so that phase corrections determined using all three properties will 
be better than those determined using only one or two of the properties. The determined 

5 properties for materials in a subject between transducer elements and a region to be 
treated are inserted into formulas developed by the inventors. The formulas use the 
determined properties and known characteristics, e.g., propagation speed, to calculate 
phase adjustments for each of the transducer elements. Compensating the phase in the 
transducer elements more properly focuses energy within the subject, permitting effective 

10 noninvasive therapy, even in sensitive areas such as the brain. Also, x-ray attenuation can 
be related to phase distortion and acoustic attenuation. Formulas relating to x-ray 
attenuation to phase distortion and acoustic attenuation can be used to compensate for 
such variances by adjusting transmitted phase and amplitude accordingly. 

Referring to FIG. 1, an ultrasound therapy system 10 includes an imager 12, a 

15 phased array 14 of n transducer elements 16, a signal adjuster 1 8, a controller 20, and a 
frequency generator 22. The system 10 is configured to determine characteristics of a 
skull 28 of a patient 30, and to apply ultrasound energy (e.g., in the range 0.01 MHz to 10 
MHz) that is focused in the patient's brain. Signals to the array are provided by a driving 
arrangement similar to that reported in Daum et ah, 'Design and Evaluation of a 

20 Feedback Based Phased Array System for Ultrasound Surgery." IEEE Trans. Ultrason. 
Ferroelectr. Freq. Control 45(2):431-4, 1998, but with a driving frequency selected 
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between about 0.1 MHz and about 10 MHz. The power and phase to each transducer 
element 16 may be manually controlled or automatically controlled using software and 
feedback. 

The array 14 of transducer elements 16 is configured to be disposed on or near the 
5 external surface of the patient's skull 28. The array 14 is configured in a curved shape 
(e.g., spherical, although sections of other shapes are possible such as planar) conducive 
for being placed on the patient's head and for focusing ultrasound energy at a distance 
from the surface of the array 14. The transducer elements 16 of the array 14 are 
piezoelectric transducer elements arranged in the array 14 as shown. The transducer 

10 elements (e.g., 1 cm 2 piezoelectric ceramic pieces) are mounted in silicone rubber or any 
other material suitable for damping the mechanical coupling between the elements 16. 
Other materials may also be used for the array construction. For example, the array 14 
may be formed from one or more pieces of piezocomposite material, or any material that 
converts electrical energy to acoustic energy. To help assure maximum power transfer to 

15 the transducer elements 16, the transducer elements 16 may be configured for electrical 
resonance at 50Q to help match input connector impedance. 

The array 14 is coupled to the signal adjuster 1 8 that is further coupled to the 
frequency generator 22. The frequency generator 22 is configured to provide a common 
radio frequency (RF) signal as the input signal to the signal adjuster 18. The radio 

20 frequency generator 22 can be of any type that will produce the appropriate signals for the 
signal adjuster 18. For example, the generator 22 may be a Model DS345 generator 
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available from Stanford Research Systems. Together, the radio frequency generator 22 
and signal adjuster 18 are configured to drive the individual transducer elements 16 of the 
array 14 at the same frequency, but at different phases (and possibly different 
amplitudes), in order to transmit ultrasound energy through the patient's skull 28 and 
5 focus the energy at a selected region within the patient's brain. The generator 22 is 
coupled to the adjuster 18 to split the generator's joutput signal to provide n input signals 
to the signal adjuster 18. 

Coupled to receive each of the n input signals from the frequency generator 22 are 
n pairs of amplifiers 24,-24 n and associated phase shifters 26,-26,, of the signal adjuster 
10 18. Each pair of phase shifter 26 and amplifier 24 represents a channel of the signal 

adjuster 18. The phase shifters 26 are configured to provide n independent output signals 
to the amplifiers 24 by altering or adjusting the phase (and possibly amplitude) of the 
incoming signals from the generator 22 by respective phase shift factors ai-a„. The 
phase shifters 26 provide approximately 1 degree precision (8-bit resolution, although 
15 lower phase resolution may be adequate for many applications). The amplifiers 24 r 24 n 
are configured to amplify the signals from the phase shifters 26 and to provide the 
amplified signals to the transducer elements 16 through connections, e.g., coaxial cables, 
individually connecting the amplifiers 24 and the transducer elements 16. An internal 
power meter is configured to monitor power supplied by the amplifiers 24. 
20 The phase shift factors cii-On of the phase shifters 26 provide steering of the 

ultrasound beam absent an object in the path of the ultrasound energy and also provide 
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compensation for phase distortion in the ultrasound output by each transducer element 16 
induced by the patient's skull. The component of each phase shift factor associated with 
steering can be computed using known techniques. The component of each phase shift 
factor (Xi-On associated with phase . distortion compensates for perturbations and 
5 distortions introduced by the skull 28, the skin/skull interface, the dura matter/skull 

interface, by variations in the skull thickness and by structural considerations such as air- 
filled or liquid-filled pockets in the skull 28. The two components that make up the phase 
shift factors aj-On are summed in order to determine the composite phase shift factors aj- 
ctn for the respective channels in order to focus ultrasound energy at a desired steering 

10 angle relative to, and distance from, the array 14. 

The phase shift factors <Xi-a„ are provided by the controller 20. The controller 20 
is logic that may be provided by software, hardware, firmware, hardwiring, or 
combinations of any of these. For example, the controller 20 can be a general purpose, or 
special purpose, digital data processor programmed with software in a conventional 

15 manner in order to provide and apply the phase shift factors aj-On to the phase shifters 26, 
although other configurations may be used. The controller 20 is configured to determine 
the phase shift factors oti-On as described below based on information obtained from the 
imager 12 as indicated by arrow 32. In particular, the information includes the thickness 
and density of the patient's skull 28 for each portion of the skull 28 between each 

20 transducer element 16 and the desired focal point in the patient's brain. Information from 
the imager 12 is conveyed directly from the imager 12 because the imager 12 is 
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configured to automatically analyze images and determine characteristics of interest from 
the images. Enough information is provided by the imager 12 to the controller 20 to 
determine the phase shift factors ai-cc,,. 

The controller 20 is configured to manipulate images from the imager 12. The 
5 controller 20 is configured to produce a 3-dimensional rendering of the patient's skull 28 
from 2-dimensional images received from the imager 12 and to determine skull thickness 
from the 3-dimensional rendering. The 3-dimensional rendering can be divided by the 
controller 20 into voxels (a volume pixel of the 3-dimensional image). 

The imager 12 is configured to obtain images of the interior of the patient's head, 

10 and in particular images that provide information regarding thickness, density, and 
structure of bone of the patient's skull 28. For example, the imager 12 may be a 
Magnetic Resonance Imaging (MRS) device or Computer Tomography (CT) device. The 
imager 12 is configured to scan the patient's skull 28 and provide information related to 
skull thickness, density and structure. This information includes 2-dimensional images of 

15 varying intensity from which 3-dimensional renderings can be made and from which 
thicknesses and densities can be determined and/or inferred. Three-dimensional image 
acquisition may also be possible and can be used. Also, if the imager 12 is a CT device, 
the imager 12 can determine and provide the CT number (also called Hounsfield number) 
for each pixel in images provided by the imager 12. 

20 Referring also to FIG. 2, the skull 28 includes two layers 50, 54 of trabecular bone 

and a layer 52 of cortical bone. For reference, a scan direction of the imager 12 is defined 
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10 



as the Cartesian x-axis so that the image plane is defined by the y-axis and z-axis. From 
the 3-D rendering, two vectors 34, 36 on the surface 38 of the skull 28 are determined 
using the 3 rd nearest neighboring points on the surface in the x and y directions from the 
point of interest ro on the skull's surface 38. The controller 20 is configured to use the 
vectors 34, 36 to calculate a vector 40 that is normal to the surface 38. The controller can 
calculate a scalar product of the vector 40 and a vector 42 that is the direction of 
propagation of energy from a transducer element 16 x . This scalar product is the incident 
angle: 



(1) 



A Cartesian (normalized) unit vector n = (u,v,w) may be calculated by the controller 20 
from this incident angle yj. Thicknesses, Dj, D2, D 3 , of each layer 50, 52, 54, 
respectively, are obtained, using intensity variations, starting at a location on the skull 
surface r 0 = (xo,yo,zo) and continuing to a point on the next layer r = (x,y,z), 

15 r= r 0 +Dn. (2) 

Near-normal propagation through homogeneous skull bone is assumed. Treating 
the layers 50, 52, 54 as a single layer for purposes of ultrasound propagation speed, the 
following simple relation can be used to express the phase shift due to propagation 
through a given point on the skull. 

20 ^360^-^) (3) 
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where f is the driving frequency, the sound speed in water is cq, the speed in the skull 
bone is c s , and D is the thickness of the skull. This convention of subtracting the phase 
shift due to direct propagation of the field in water from the shift due to the skull 
produces a positive phase shift 
5 The three layers 50, 52, 54 can be treated separately to determine phase 

compensation due to the thicknesses D b D2, D 3 of the layers 50, 52, 54 as part of a three- 
layer model. For the three-layer model, it has been assumed that the skull 28 consists of 
individual homogeneous layers. The speed of sound is assumed to be 2.5 x 10 3 m/s for 
the central layer 52 and 2.9 x 10 3 m/s for the inner and outer layers 50, 54. The expected 
10 phase shift across the skull 28 using this three-layer model is: 

# = 360fiD n (±-±) (4) 

with c„ equal to the speed of sound in the n^ layer and D n being the thickness of the n* 
layer. 

It has been found that density of the skull 28 affects the propagation through the 
15 skull 28. Consequently, the phase of ultrasound signals passing through the skull 28 is 
affected by the skull's density, and phase correction based on the skull's density can help 
improve focusing of ultrasound energy from the array 14. 

To assess the contribution of internal density variation to the deviation from the 
expected single-layer homogeneous values calculated from Eq. (3), a correlation between 
20 the mean CT image intensity along the ultrasound axis of propagation 42 and 

experimental error is used. Using a 3-d image constructed from 2-d CT images, each 
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0.15 mm 2 voxel is assigned an intensity value. It is assumed that the intensity is linearly 
proportional to bone density and density is scaled to MKS units using air and water in the 
image as reference intensities. Mean intensity is determined by summing the CT intensity 
values along the axis of propagation 42 inside the bone and dividing by the total number 
5 summed of voxels. The voxels may include air-filled or liquid-filled pockets. The sound 
speed for such voxels is assumed to be the speed of sound in water, or if air-filled, then 
complete reflection of the ultrasound can be assumed. 

Error due to skull density has been calculated as the difference between the 
measured phase and that given by Eq. (3) (for a single-layer model). An empirical 
10 correction factor has been obtained by fitting (using a polynomial curve fit) percent error 
as a function of the mean intensity. The correction factor is: 

C{p)=f J A n p N -'' (5) 

where there are N+l coefficients A n to be determined for an N m degree fit, and p is the 
density. Applying a 3 rd order fit using a 1 17 point data set and five skulls at 0.51 MHz 
15 yielded coefficients, in MKS units, of: 
Table 1 



Ao 


1.1424*408 


A, 


-7.5377e-005 


A 2 


0.1645 


A 3 


-118.689 
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These coefficients may be used in Eq. (5), and Eq. (5) applied to Eq. (3) as part of the 
single-layer homogeneous model to yield a corrected phase shift of: 

(6 > 

where <t> contains the initial phase values obtained using Eq (3). 
5 The density and thickness measurements may also be applied toward adjustment 

of the skull sound speed. To adjust the speed value for a single layer, the controller 20 is 
configured to fit the speed as a function of density according to: 

where <f>(p) is the measured phase shift as a function of density. A polynomial fit as 
10 shown in Eq. (5) is used by the controller 20 to find the speed of sound values. For the 
three layer model, two sound speeds are calculated. These two speeds are the speed cj for 
the trabecular layers and the speed c a of the cortical (central) bone. Given the thicknesses 
D b D 2 , and D 3 , respectively, for the three bone layers 50, 52, 54, the speed c s of the 
cortical layer is fit by the controller 20 as a function of density according to: 

c 0 c u 360JD 

The polynomial fit for c f (p) is performed using Eq. (5) over a series of trial functions for 

c iS . The final sound speeds are the c f (p) and c u (p) that most closely correlate by 

standard deviation or other methods with direct skull measurements obtained by fitting a 
large skull sample. 
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Referring to FIGS. 1 and 3, in operation, the controller 20 uses information from 
the imager 12 to adjust the phase of the transducer elements 16. At stage 60, the imager 
'12 takes images of the patient's skull 28. This may be done remotely in space, e.g., in 
another hospital, from the controller 20. At stage 62, the imager transmits information, 
5 e.g., intensities of portions of the image, related to skull thickness, density, and structure 
to the controller 20. This transmission can be separated in time from when the imager 12 
takes the image and can also be performed by human intervention, e.g., by recording the 
images on a CD and replaying the CD in the controller 20. At stage 64, the controller 20 
manipulates this information to determine thicknesses of layers, and to identify the layers, 

10 of the patient's skull 28 by analyzing intensities of portions of the image. At stage 66, the 
controller 20 determines excitation correction factors such as the phase shift factors otpOn 
(or amplitude corrections as discussed below) and transmits the phase shift factors a r On 
to the adjuster 1 8. At stage 68, the frequency generator supplies energy to the adjuster 18. 
The adjuster 18 adjusts the phase of the energy from the frequency generator 22 

15 according to the phase shift factors ai-cin from the controller 20. The adjusted energies 
are sent to the transducer elements 16 that convert the energies into ultrasound waves and 
transmit these waves into the patient's skull 28. Images from the imager 12 can be 
supplied in real time and the amplitude and/or phasing of the energy applied to the patient 
30 changed in response to the images. 

20 The power supplied to the patient's skull 28 depends on the type of therapy. For 

ablation, approximately 2-3 kW for approximately 1 0 seconds using 64 transducer 



25 



WO 03/01 7843 PCT/US02/27130 



elements may be used. If more transducer elements 16 are used, then less total power 
may be used and vice versa. For opening the blood-brain barrier, about 100 times less 
power than for ablation may be used due to preformed gas bubbles in the area of interest. 
The ablation power can also be reduced by the preformed gas bubbles. Using bursts of 
5 energy has been found to reduce, if not eliminate, affects on phase due to standing waves 
that may occur if the transducer elements 16 constantly emit energy. 

Other embodiments are within the scope and spirit of the appended claims. For 
example, the array 14 of transducer elements 16 shown in FIG. 1 may contain fewer 
transducer elements 16 than shown. The phase shift factors apotn may be pre-stored in 

10 the channels of signal adjuster 18 instead of being provided by the controller 20. Further, 
functions described as being performed by the controller 20 could be performed by the 
imager 12 and vice versa, or by a person using the system 10, e.g., calculating densities 
and providing input data to the controller 20 regarding phase shift factors a r «n. 

Referring to FIGS. 4-5, a noninvasive process 70 for obtaining phase shift factors 

15 oti-On using a system 90 operates the transducer elements 1 6 in both transmit and receive 
modes. In the system 90, a controller 92 is similar to the controller 20 (FIG. 1), but is 
also configured to control switches 94 r 94 n such that the transducer elements 16 are 
connected to the controller 92 in a receive mode and to the adjuster 1 8 in a transmit 
mode. 

20 At stage 72, the controller 92 controls the switches 94 to be positioned such that 

one switch, e.g., switch 94i, is in the transmit mode. The frequency generator 22 
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produces a low-frequency ultrasound wave train of about 5-30 cycles. Some of this wave 
train reflects off of the skull 28 and some of the wave train passes through the skull 28 
into the patient's brain. 

At stage 74, the controller 92 causes the switch 94\ to switch into the receive 
5 mode, connecting the transducer element 16\ to the controller 92. The wave train reflects 
off of the skull 28 and reflections are received by the transducer element 16\ and recorded 
by the controller 92. The propagation speeds and the pulse train lengths are such that 
energy from a pulse train is still being transmitted into, and reflected by, the outer surface 
of the skull 28 when energy from that pulse train that was reflected by the inner surface of 
10 the skull 28 passes through the skull 28 toward the array 14. Thus, the received 
reflections of the inner and outer surfaces are superimposed. 

At stage 76, a query is made as to whether a series of frequencies to be transmitted 
has been completed. If not, then the process 70 returns to stage 72 to transmit a different 
frequency in the series. If the series has been completed, with corresponding recorded 
15 amplitude data as a function of frequency, then the process 70 proceeds to stage 78. 

At stage 78, the controller 92 analyzes and processes the recorded information to 
determine the phase shift factors otpcc*. The controller 92 deconvolves the recorded 
frequency response that includes the superimposed reflections from the interfaces at the 
inner and outer surfaces of the skull. The superimposed reflections produce a periodic 
20 appearance of local maxima and minima in the data as a function of frequency. 
Information from images, e.g., CT or MRI images, are used by the controller 92 to 
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determine the locations of the inner and outer surfaces, and thus the skull thickness. The 
controller 92 estimates the phase shift (A<|>) using the distance, Af, between like extrema, 
e.g., peaks, according to: 

A<|> = 27tf [ (d/c 0 ) - (l/(2Af)) ] (9) 
5 where f is a frequency to be transmitted for which the phase correction is to be used, d is 
the skull thickness determined from the image, and Co is the speed of sound in water. 
Adjustments may be made to this phase shift using density correction techniques. The 
process 70 is repeated for the remainder of the transducer elements 16. 

Phase shifts and acoustic attenuation due to transskull sonication can also be 

10 determined from analysis of CT images. For purposes of determining phase shift and 
acoustic attenuation, it can be assumed that there is essentially a one-to-one correlation 
between x-ray attenuation coefficient |i and the speed of sound c in bone and between \x 
and acoustic attenuation. 

Equations relating x-ray attenuation to phase and acoustic attenuation have been 

15 developed as discussed below, with the equations including parameters determined from 
fitting to experimental data. Experimental data has been determined by applying 
ultrasound to a pig skull sample using a 1" transducer element made by Panametrics of 
Waltham, Massachusetts, operating at 730 kHz and positioned about 7 cm from a 0.6mm 
hydrophone. The skull was attached to a 3-D positioning system and the transducer 

20 element was driven by a 1 0-cycle sinusoid signal produced by a waveform generator such 
as a Wavetek model 395 made by Wavetek Ltd. of Norwich, United Kingdom, or a 
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waveform generator made by, e.g., Fluke Corporation of Everett, Washington. The 
hydrophone response was recorded on an HP54602 digital oscilloscope made by Hewlett 
Packard of Palo Alto, CA. Measurements of the skull sample were obtained by 
immersing the sample in water and using dual-slice spiral acquisition with 0.5 mm slices. 
5 Phase shift through the skull can be related to x-ray attenuation by assuming a 

formulaic relationship and fitting actual data to the formula. The x-ray attenuation 
coefficient \x is related to the CT number CT n according to: 

^ = k(CT n +1,000) (10) 
where k is a well-known physical constant. The units of \i may be chosen such that k=l . 

10 The phase shift due to the passage through the bone is given by: 
A(p=27rff(l/c(x)-l/c 0 )dx (1 1) 

where f is the frequency, c(x) is the speed of sound at point x and Co the speed of sound in 
water (or soft tissue, that is practically the same). The integration in Eq. 1 1 is performed 
along the line of propagation of the ultrasound piergy. The functional dependence on \i is 

15 assumed to be: 

l/c(x) = A' + B'(l/K*)) + C'OMx)) 2 (12) 
where A', B' and C are parameters fitted to the experimental data determined from 
experiments using the setup described above. Approximations other than that shown in 
Eq. 12 may be used. FromEqs. 1 1 and 12: 
20 A(p = A(p 2 - Pl +l) + BS(l/^(x)) + Ci:(l^(x)) 2 (13) 

where A = 27tf(A' - l/c 0 ), B = 27ifB' , C = 27cfC', and p2, pi are the first and last pixels 
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within the bone along the line of propagation and the sums are performed within these 
boundaries. The values of P2, pi are obtained manually from intensity profiles drawn 
along the lines of propagation. These values may also be obtained automatically using 
well-known edge-detection techniques. 
5 The equation for A<p can be modified to include a thickness-independent 

parameter. Thus, Eq. 13 can be modified to: 

Aq > = A(p2-pi+l) + BS (l/n(x)) + C2 (l/^x)) 2 + J (14) 
where J is the phase shift at zero thickness and its non-zero value is probably due to 
experimental error. The value of A9 can be used to adjust the phases of elements of an 
10 array of transducer elements such as the array 14 of the system 10 shown in FIG. 1. 

From experiments using the setup described above, exemplary values of the 
parameters A, B, C, and J were obtained. These values were A= 0. 1 126, B= -0.0610, C= 
-0.0759, and J= -0.1375. Using these values in Eq. 14 yielded superior results (e.g., a 
standard deviation of 21°) when compared to A(p = 2rcfx (1/ c b -1/co) and assuming c b = 
15 2650 m/sec and Co = 1 500 m/sec (e.g., a standard deviation of 54°). 

Acoustic attenuation through the skull can also be related to x-ray attenuation by 
assuming a formulaic relationship and fitting actual data to the formula. The acoustic 
attenuation through bone is given by: 

ln(W 0 ) = - W;(x)dx (15) 
20 Where I and I 0 are the acoustic intensities before and after the passage through the bone 
respectively and |i AC (x) is the local acoustic attenuation coefficient. The relation 
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between the acoustic and x-ray attenuation coefficients can be approximated as: 

*i AC (x) = E+F m(x) +G(fx(x)) 2 (16) 
Approximations other than that shown in Eq. 16 may also be used. Integrating over the 
set of pixels along the line of propagation, as explained above, yields: 
5 -ln(I/I 0 ) = E(p2-pi+l) + FE^ (x) + G2 (tfx)) 2 + H (17) 

where H is a thickness-independent parameter thaf reflects the fraction of the intensity 
reflected at the interfaces (air-to-skull, between skull layers, and skull-to-brain). Again, 
from experimental results using the setup discussed above, acceptable values for the 
parameters E, F, G, and H were determined to be E= 0.7552, F= -0.9007, G= 0.2946, and 
10 H= -0.9 121. Using information about the acoustic attenuation, the magnitude of 
transmitted energy can be adjusted using the amplifiers 24 based on signals from the 
controller 20. 

In operation, the imager 12 supplies the CT numbers of images to the controller 
20 that uses the CT numbers and the experimentally-determined values for A, B, C, J, E, 
15 F, G, & H to determine the phase correction factors ai-<Xn and amplitude adjustment 
factors for the amplifiers 24. The system otherwise operates generally as described 
above, with the adjuster 18 adjusting phase and amplitude of energy from the generator 
22 as indicated by the controller 20. 

Still further embodiments are within the scope and spirit of the appended claims. 
20 For example, while the above discussion focused on ultrasound energy transmitted 

through human skulls, other frequencies of energy and other subjects may be used, and 
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other areas of subjects may be targeted for receiving energy. For example, non-human 
subjects are acceptable, as well as transmitting energy into areas of human or non-human 
subjects aside from through the skull and into the brain. Also, other techniques for 
determining and/or controlling the phase and amplitude of signals transmitted into the 
5 patient 30 may be used. For example, techniques that may be used include those of 
concurrently-filed U.S. Patent Application Ser. No. [Not Yet Assigned, filed November 
28, 2000; Attorney Docket Number 254/187] entitled "SYSTEMS AND METHODS 
FOR FOCUSSING AN ACOUSTIC ENERGY BEAM TRANSMITTED THROUGH A 
SKULL," and concurrently-filed U.S. Patent Application Ser. No. [Not Yet Assigned, 

10 filed November 28, 2000; Attorney Docket Number 254/189] entitled "SYSTEMS AND 
METHODS FOR CONTROLLING A PHASED ARRAY FOCUSED ULTRASOUND 
SYSTEM," both of which are incorporated here by reference. 

Further embodiments are also within the scope and spirit of the invention. For 
example, the system 10 shown in FIG. 1 can be configured to apply a projection 

15 algorithm to identify phase and/or amplitude adjustments to be applied by the amplifiers 
24 and/or the phase shifers 26. The algorithm is based on a layered wavevector- 
frequency domain model that propagates ultrasound from a hemispherical transducer 
array through the skull 28 and uses input from CT scans of the skull 28. The projection 
algorithm, that the controller 20 is configured to perform, calculates driving phases and/or 

20 amplitudes for each transducer 16 in the array 14 to improve, an possibly maximize, the 
signal at an intended focus. Experimental results have shown that this technique can be 
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used for completely non-invasive ultrasound brain therapy, including surgery. 

The projection algorithm as implemented by the system 10, and especially the 
controller 20, concentrates on a noninvasive focusing method that uses CT images to 
predict how the ultrasound will propagate through the skull 28. The algorithm predicts 
5 the behavior of the ultrasound field after passing through the skull bone 28, which causes 
significant reflection, diffraction and absorption of -the field. To implement the 
algorithm, the system 10 accurately obtains knowledge of the thickness and internal 
structure of the skull bone 28 and achieves precise registration between the skull 28 and 
the ultrasound array 14. The algorithm provides a relatively uncomplicated, 

10 computationally-feasible mechanism for determining phase and/or amplitude corrections 
without oversimplification. 

The projection algorithm uses thickness, density, and orientation information 
obtained from CT images of the head 28. This information is entered into a propagation 
model that can propagate ultrasound through arbitrarily oriented (non parallel) layers. 

15 Operating in the wavevector-frequency domain, this model rapidly propagates sections of 
the ultrasound field through the skull 28. Two specialized cases of the model are 
examined for the present focusing algorithm. The first case propagates ultrasound 
through the skull 28 with normal incidence, and assuming that the skull is homogeneous. 
The second case considers the angle between the skull surface and the incident acoustic 

20 field. Both of these models have been tested using ex vivo human skulls and hydrophone 
measurements in a water bath. 
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The Projection Algorithm 

A projection algorithm in the wavevector-frequency domain is used to determine 
propagation through the skull 28 and correction factors for the ultrasound array 
5 transducers 16. The projection algorithm assumes that the acoustic field satisfies the 
linearized Stake's equation, 

[l+r|]v 2 j p(r,0 = p-^0, (18) 

where x is the relaxation constant for the medium and c is the real sound speed. While 
therapeutic arrays operate at high powers, the linear approximation is found to be 
10 appropriate for propagation outside of the focus. A Helmholtz equation in Cartesian 
coordinates is obtained by the substitution of the Fourier integral 

PirA-^fiMk^z^e^dk^da, (19) 

into Eq (1 8), where the Cartesian wavenumbers are given by k x and k y and (pis the angular 
frequency. The resulting equation, 

15 (V 2 +P) p(k xy k y ,z,a>) = 0 , (20) 
is expressed in terms of a complex wavenumber k given by 



i=J4r-V— V-V- < 21 > 

V c l+ira> y 

To propagate the field in a plane at a distance zo in front of a source to a new plane z, the 
advanced solution of Eq. (20) is used: 
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(22) 



10 



15 



20 



In wavevector space, the field recorded in a plane zo is thus related to the field at any 
other plane z by a simple transfer function given in the right hand side of Eq. (22). The 
pressure field in space at z may b t e obtained by taking the inverse Fourier transform of the 
equation. 

Propagation in wavevector space can be performed using Eq. (22) after the spatial 
Fourier transform of the field is calculated for a given plane. This planar cross section of 
the field, zo, is chosen to be planar and perpendicular to the propagation axis. The field 
can be projected to an arbitrary new plane z. The array 14 is sectioned into small areas 
that are individually propagated through the skull 28 in a manner that, as described 
below, yields more easily in Cartesian space. The controller 20 divides the array 14 into 
segments equal to the array's transducer elements (e.g., 280 mm 2 ), although this does not 
have to be the case. For example, a large element could be segmented into several 
smaller elements. 

The controller 20 calculates the acoustic pressure in water due to each segment 
over a plane located directly above the skull's outer surface. Using the backward 
projection approach, the controller 20 first calculates the field of each transducer element 
16 at the transducer's geometric focus in water. The field at the focus is calculated using 
the Rayleigh-Somerfeld integral, 
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over a planar area 20mm x 20mm perpendicular to the geometric center of the transducer 
16. A water sound speed of c = 1488 m/s and density of p = 1000 kg/m 3 is used along 
with a normalized surface displacement velocity v . The controller 20 preferably 
calculates this field only once for a particular transducer, and uses it repeatedly. 
5 Alternatively, the effects of imperfections on the array could be added by replacing Eq. 
(23) with physical measurements of the field over this area for each element 16 in the 
array 14. 

The field due to the n-th element is Fourier transformed by the controller 20 
giving a planar field p{k x , k y9 a>,z 0 ) . This planar field is rotated in the wavevector- 
10 frequency domain, defining a new projection plane, z 0 ■ . Referring to FIG. 6, the spatial 
Fourier transform method, sometimes referred to as the angular spectrum approach, may 
be viewed as an angular decomposition of the field. That is, each point describes a planar 
wave with frequency co passing through the measurement plane at an angle equal to 

a x = sin" 1 k x — from the z-axis in the x-z plane and a v = sin" 1 k v — in the y-z plane. The 

a> y y a> 

15 amplitude of the wave is equal to the modulus of p . FIG. 6 illustrates this concept for a 
two-dimensional case. In three dimensions, the propagation angles may be expressed in 
terms of a polar angle 9 and azimuthal angle <p: 
tana x = tantfcos^ 



tana^ = tan (9 sin ^ 



(24) 



For a transient signal, each temporal frequency possesses a unique angular spectrum, 
20 forming a three-dimensional space consisting of two spatial frequency dimensions and 
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one temporal frequency dimension. The tested signals are harmonic and may be 
expressed in two dimensions. 

Viewing the wavevector plane as an angular decomposition, to the controller 20 
can perform a passive coordinate rotation. To orient the plane along a new projection 
5 axis, z' 0 , each point p in wavevector space is mapped to the rotated plane such that: 



where a x is the rotation angle in the x-z plane and a y is the rotation angle in the y-z 

plane. These angles are be obtained by inverting the polar and azimuthal orientation of 
the element using Eq. (24). Projection along the new axis to an arbitrary new plane z' is 
10 achieved using Eq. (22). The inverse Fourier transform performed after this mapping 
provides the pressure field in the plane rotated to {0*°* ^ R0T ^ , where the spherical polar 

angles are related to a x ROT and by Eq. (24). Transskull propagation begins by 
projecting to the rotated plane at z' 0 to the plane z' directly above the skull's outer surface. 



15 not known a priori, fee projection approach will yield their locations. The controller 20 
obtains the angles directly from the Fourier transform of the pressure field, taking angles 
associated with the maximum amplitude of p(k x9 k yJ (D 9 z 0 ) as the propagation axis. 




(25) 




For the physical array 14, where the exact orientation {0 7 <f>) of the elements 16 is 



Propagation through layered media in wavevector space can be achieved by 



37 



WO 03/017843 



PCT/US02/27130 



determining exactly how the layers distort each point in the space. For an arbitrary field 
incident upon a series of planar layers, each wavevector component over the plane at zq is 
associated with a wavevector k 0xy (tt?>=k x (fl?)i+k y (a?)j , where the indices (x,y) indicate 
the position of the component in k-space. The angle of incidence between a given 
5 wavevector and the vector normal to the outer surface of the n* layer is given by {&) , 

and the transmitted angle by y tmf {®) . For a layer with sound speeds c n and density p n , 
the pressure transmission coefficient at its lower interface is given by 

l my\®) p^COS^+p^C^COS^ ' ( 26 ) 

where the incident and transmitted angles yj, y t are understood to be a function of 
10 frequency. Since the problem has been decomposed into a series of plane waves, 

transmission of the field from the outer surface to the next layer may be treated either by 
using a ray method, or by applying the continuity of the pressure and normal particle 
velocity at the interface. The ray method is described below in detail, as it leads to the 
general case of propagation through non-parallel layers. In this first-order propagation, 
15 multiple reflections within the layers are neglected, although such higher order reflections 
could readily be added. 

Regardless of whether or not the layers are parallel, the transfer function may 
readily be written in a closed form. In addition to the thickness across the z-axis, Zn, the 
sound speed c n , and density of each layer Pn, the unit vectors normal to the layer surfaces 
20 h n are calculated. For a given initial wavevector k^ , the ray path from (0,0, zo) between 
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any two surface interfaces traverses a distance of 

pq. fc.-o-^ (27) 

******* 

where, as depicted in FIG. 7, r^is the vector extending along the layer from z-axis to the 
intercept of the layer with the ray. The unit vector along the wavevector's path is given 
5 by . Again, the frequency dependence on the wavevector orientation is understood. 

It follows that the ray position vector is = |^|^ . Although the initial wavevector 

orientation, k Qxy , is known, the direction of the wave vector in the first and subsequent 
layers is calculated using the relation 

— (" n+ i x ) = — (>L * Lb, ) > ( 28 > 

10 which is a consequence of Snell's law in three dimensional space in which the incident 
wavevector, the transmitted wavevector and the normal vector all be in the same plane. 
The transmitted wavevector on the right hand side of Eq. (28) may be obtained by 
crossing both sides of the equation with h n+i . Using cross product relations, it may be 
shown that 

15 Llxy=—{kr*y -COS^^, J + COS^H^ (29) 
Si 

The incident unit wavevector of the n* layer is equal to the transmitted wave of the (n-l)* 
layer. With the exception of the 0 th layer shown in FIG. 7, is also found in order to 

calculate Eq. (27). Given the thickness across the z-axis of each layer, the point of 
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intersection of ray R n with the surface of the nH layer is 

7 n+\=K+r n -z»' (30) 

Over a series of N layers, the phase of a ray reaching the plane ^(k^co) is 

the sum of the phase contributions over each path length given by Eq. (27). The spatial 
5 phase at the plane z, is related to the ray phase at N by 

&(*j^^) = &^ 7. A ray leaving the 

initial plane with a polar angle y 0xy (a) will arrive at the plane z with a new orientation 

Ynv {*>) determined by the N-l unit vector k N _ t . The spatial phase at plane z is related 
to the spatial phase at plane z Q by 

10 tfo,.*.*)-^,*.*,^*. ft-*)-*" -sinr^H, (31) 



given 0(k Ov ,a>,z o ) , the phase of p(k x ,k y ,a>,z 0 )at the initial plane. 

The pressure over the plane at z can be expressed in terms of the ray phase 
presented in Eq. (31) and the transmission coefficient in Eq. (26) by 

Kk^z) = p(^^,z 0 K'"^ W ^n^(«>* ( ^- ffl,z) j , (32) 

15 where the terms in square brackets on the left hand side of Eq. (32) can be viewed as an 
operator that maps the field from ko-space to a new k N -space. In practice, this mapping 
uses interpolation to produce a linearly spaced matrix at z. 
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Non-Planar Surfaces 

The layered-projection method may be used for the propagation of a field through 
a curved surface, provided the surface is sufficiently smooth relative to the highest 
relevant wavenumber. The field is projected to a plane near the surface and is divided 
5 into a series of virtual sources. The region of the surface penetrated by a given source is 
approximated as planar, with the actual surface over the beamwidth of the source 
preferably agreeing with this approximation to at least within V4. of the maximum 
frequency k m . At the same time the source diameter S is preferably large enough that 
the smaller frequency components of the source are in the near field relative to the 



10 Rayleigh distance, with S > 




, where is the smallest appreciable frequency of 



the signal and z is the distance from the starting plane to the surface. Although the planar 

field is uniquely divided into a continuous series of sources, the fields of the sources will 

overlap on the surface. 

When propagating through multiple layers, each of the surfaces is segmented as 
15 described above. It is possible that components with high transmission angles will leave 

the surface of interest. 

Thus, multiple layers of the skull 28 can be taken into account by determining the 

transmission angle through each layer, using this angle and layer thickness and orientation 

(and the next layer's orientation) to determine length of propagation through the layer, 
20 and applying the density to determine phase and attenuation (as discussed below) through 

each layer based on the propagation length through each layer. 
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Registration of Images and Array 

Standard techniques may be used to register the array 14 and images of the skull 
28. Registering the array 14 and the images helps to ensure that the proper focal point for 
5 the ultrasound is achieved. Proper registration helps ensure that the proper portion of the 
skull 28 is being analyzed to yield the accurate thicknesses, densities, orientations (e.g., 
incident angles), and locations of layers of the skull 28 between the array 14 and the 
desired focal point within the skull 28. With accurate measurements, the controller 20 
can determine accurate phase and amplitude corrections for the amplifiers 24 and the 

10 phase shifters 26 to focus ultrasound from the array 14 at a desired location inside the 
skull 28. Standard techniques for registering the array 14 and the skull images include 
using a reference frame that attaches to the array 14 and to the skull 28 and that has 
reference markers that can be identified in CT and MR images. The reference markers 
can be used to rotate CT and MR images to register the images (with CT images yielding 

15 bone characteristics, while MR images may be used during therapy due to better soft- 
tissue resolution). Alternatively, one or both of the CT and MR images may be rotated 
until they match. Other techniques that yield registration so that the array's orientation 
relative to the skull 28 can be determined are also acceptable. 

20 Skull Sound Speed 

Based on a study of over 1000 measurements performed with 10 skulls using a 
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0.51 MHz transducer, a linear fit of sound speed as a function of skull density was 
obtained. A first-order (linear) fit obtained using Eq. (22) yields the equation: 

c cff =2.06^p-1.54x 10 3 f (33) 

in MKS units over the range of thickness-averaged densities between 1.82 x 10 3 kg/m 3 
5 and 2.45 xlO 3 kg/m 3 . Density measurements obtained by analyzing the CT images are 
input to Eq. (33) to provide the effective skull sound speed values in the focusing 
algorithm. Eq. (33) is exemplary only, and not limiting. Eq. (33) is based on 
experimental data, and thus dependent on the acquired data. Further, other relationships 
of sound speed and density are possible, and acceptable, including other linear data fits. 
10 Also, other fits, e.g., non-linear fits, of sound speed as a function of skull density may be 
used. 

Skull Attenuation 

Experimental data can be, and has been, used to determine a relationship between 
15 ultrasound attenuation and skull density. It was found that attenuation increases as the 
frequency increases and generally decreases nonlinearly as the density increases. Further, 
the relationship of attenuation and density has a minimum and a slightly increasing 
attenuation when the density increases beyond approximately 2000 kg/m 3 . 

20 Focusing Algorithm 

To produce a focus through the skull 28, the acoustic pressure is calculated by the 
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controller 20 at the intended focal point for each element 16 after propagation through the 
skull 28. The controller 20 can perform its calculations under at least two different 
scenarios. In the first case (Case 1), the controller 20 ignores incident angles and assumes 
the ultrasound is normally incident upon the skull 28. In the second case (Case 2), the 
5 controller 20 uses the incident angle y\ determined from Eq. (26) over the skull surface to 
include incidence in the model. The pressure phase is compared with the phase expected 
if the skull 28 were not present The change in phase caused by the skull 28 is recorded 
and used for correcting the driving phase of the transducer array 14. The driving phase of 
each element 16 is adjusted by an amount 
10 W = *rg(P(r)/P 0 (r)), (34) 

where P is the acoustic pressure at the point in the skull 28 and P 0 is the acoustic pressure 
expected at the same point if the skull 28 were not present. Similarly, a large element 16 
could be divided into M sections for propagation, with the driving phase is adjusted by 

^ = arg|fp B (r)/? 0 (r)j. (35) 

15 

Operation of the system 10 implementing the spectral wavevector algorithm 
proceeds similarly to the process shown in FIG. 3. The imager 12 takes a CT image of 
the skull 28. From this image, the controller 20 determines locations, orientations, 
thicknesses, and densities of layer portions between respective array elements 16 and a 
20 desired focal point inside the skull 28. The controller 20 can assume that the skull 28 is a 
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single layer or can take into account different skull layers. From the determined skull 
data, the controller 20 determines lengths traveled by ultrasound from the elements 16 
through the skull 28 toward the focal point. The controller 20 applies relationships of 
phase and amplitude to density, either taking into account incident angles or assuming 
5 planar propagation. The controller 20 thus determines propagation speed and attenuation 
(e.g., using Eq. (33) for propagation speed and an experimentally-derived relationship for 
attenuation) through the skull 28, and from these data and length traveled through the 
skull 28 determines the phase and amplitude deltas for each array element 16. With these 
deltas, the controller 20 determines phase and amplitude correction factors to be applied 
10 by the phase shifters 26 and the amplifiers 24 to yield focused ultrasound at a desired 
location inside the skull 28. 



Experimental Data 

The entire projection algorithm was implemented in Matlab, on a series of 
15 Pentium and AMD-based computers over an NT network. A typical processing time for 
producing a focus within the skull 28 for 320 elements 16 using a 1 GHz machine with 
256 MB of RAM was approximately five hours. 

The focusing method was tested on ten Formalin fixed ex vivo human skulls. 
Acoustic properties of the fixed skulls are assumed to be similar to fresh skulls. CT 
20 Images of the skulls were obtained and used in the focusing algorithm, treating the skull 
bone as homogeneous, so that in Eq. (31) N=3. Although the wavevector method is 
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capable of compensating for both the magnitude and phase of the driving signal, only 
phase correction was considered. The transducer array was a 30-cm-diameter, 0.74 MHz 
hemisphere divided into 500 individual elements. Only the lower 320 elements of the 
array were used for the measurements in order to allow a larger range in the movement of 
5 the skull. A reference frame was attached to the skull to assist with registration of the 
array and the skull. A hydrophone was placed inside the skull to measure the ultrasound 
inside the skull from the array. 

For the 1 0 skulls, the peak signal without phase correction was found to occur at a 
mean distance of 1 . 1 mm from the intended focal location with values ranging from 0.45 

10 mm to 2.5 mm and a standard deviation (STD) of 0.64 mm. In 7 of the 10 cases studied, 
the main focus was scattered into two or more foci, so the peak signal refers to that of the 
largest (larger) focus. When the phasing algorithm was used, a single focus was always 
produced, regardless of whether the incident angle was considered. When the phasing 
algorithm was applied without correction for incident angle, the mean distance from the 

15 peak was reduced to 0.59 mm with values ranging from 0.22 mm to 1 . 1 mm and a STD of 
0.26 mm. With the angle of incidence included in the algorithm, the mean distance was 
further reduced to 0.48 mm over a range of distances from 0.20 mm to 0.82 mm and a 
STD of 0.22 mm. 

Peak voltage-response-squared measurements were also observed to improve with 
20 the phasing algorithm. Without phasing, the average peak was found to be 34 % of the 
value obtained by the invasive hydrophone-phasing technique with a variation from 14 % 
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to 58 % over the samples. The lower amplitudes corresponded to skulls that scattered the 
field into multiple foci. The phasing algorithm, without correction for incident angle, 
improved the peak value on average to 46% of the hydrophone-phased signal with a range 
between 22% and 58%. When the incident angle was included, the average value was 
5 45% with a range from 29% to 59%. The peak values between the two models were, on 
average, similar. However, the variance was reduced when incident angles were 
included. 

Both models significantly improved the ratio between the peak amplitude and the 
largest secondary lobes. To assure all secondary lobes were taken into account, values 
10 were obtained from searching within the 5 12 mm 3 volume obtained using the harmonic 
projection technique. Without phase correction the secondary peak amplitude averaged 
63% of the primary peak. Using the correction algorithm, secondary peaks averaged 29% 
without angle correction and 32% with the angle included. 

Still further embodiments are within the scope and spirit of the invention and the 
15 appended claims. 

What is claimed is: 
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CLAIMS 

1 . A method of delivering ultrasound signals, the method comprising: 
providing an image of at least a portion of a subject intended to receive ultrasound 

5 signals between sources of the ultrasound signals and a desired region of the subject for 
receiving focused ultrasound signals; 

identifying, from the image, physical characteristics of portions of material 
between the sources and the desired region; and 

using a spectral propagation analysis to determine at least one of phase corrections 
10 and amplitude corrections for the sources depending on the respective physical 

characteristics of the portions of the material disposed between each source and the 
desired region. 

2. The method of claim 1 wherein the physical characteristics include density 
15 and thickness. 

3. The method of claim 1 wherein the identifying identifies the physical 
characteristics of different layers of the portions of the material. 

20 4. The method of claim 3 wherein the physical characteristics include angle 

of a layer interface relative to a reference. 
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5. The method of claim 1 wherein the spectral propagation analysis includes 
a wavevector-frequency domain projection algorithm. 

6. The method of claim 5 wherein the using includes comparing projected 

5 phases from the sources to the desired region accounting for the material and ignoring the 
material. 

7. The method of claim 6 wherein the using assumes perpendicular incidence 
of ultrasound. 

10 

8. The method of claim 6 further comprising determining an incident angle 
of ultrasound to the material. 

9. The method of claim 8 wherein the identifying identifies the physical 
15 characteristics of different layers of the portions of the material and the determining an 

incident angle determines incident angles for each layer of the material. 

10. A system for delivering ultrasound signals, the system comprising: 
an apparatus configured to analyze an image of at least a portion of a subject 

20 intended to receive ultrasound signals between sources of the ultrasound signals and a 
desired region of the subject for receiving focused ultrasound signals, the apparatus 
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configured to determine, from the image using a spectral projection analysis, information 
about the at least a portion of the subject; and 

an array of sources of ultrasound signals having at least respective ones of their 
relative phases and their amplitudes set in accordance with the information about the at 
5 least a portion of the subject provided by the apparatus. 

1 1 . The system of claim 10 wherein the apparatus comprises means for 
identifying physical characteristics of multiple layers of the portion of the subject, the 
physical characteristics including density and thickness. 

10 

12. The system of claim 1 1 wherein the physical characteristics further include 

i 

incident angles of the layers relative to respective sources, the apparatus further 
comprising means for comparing projected phases of the signals from the sources to the 
desired region accounting for the portion of the subject and ignoring the portion of the 
15 subject. 

1 3 . The system of claim 1 0 wherein the apparatus comprises means for 
performing a wavevector-frequency domain projection algorithm. 
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14. A computer program product residing on a computer readable medium and 
comprising computer-readable, computer executable instructions for causing a computer 
to: 

analyze an image of at least a portion of a subject to receive ultrasound signals 
5 between sources of the ultrasound signals and a desired region of the subject for receiving 
focused ultrasound signals to identify, from the image, physical characteristics of 
respective portions of material between the sources and the desired region; and 

determine, using a wavevector-frequency domain projection algorithm, at least 
one of phase corrections and amplitude corrections for the sources depending on 
10 respective physical characteristics of portions of the material disposed between each 
source and the desired region. 

15. The computer program product of claim 14 wherein the physical 
characteristics include density and thickness. 

15 

16. The computer program product of claim 14 wherein the instructions for 
causing the computer to analyze cause the computer to identify the physical 
characteristics of different layers of the portions of the material. 

20 1 7. The computer program product of claim 1 6 wherein the physical 

characteristics include angles of layer interfaces relative to the sources. 
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1 8. The computer program product of claim 14 wherein the instructions for 
causing the computer to analyze cause the computer to compare projected phases from 
the sources to the desired region accounting for the portions of the material and ignoring 

5 the portions of the material. 

19. The computer program product of claim 1 8 further comprising instructions 
for causing the computer to determine the projected phases accounting for the portions of 
the material assuming perpendicular incidence of ultrasound. 

10 

20. The computer program product of claim 1 8 further comprising instructions 
for causing the computer to determine the projected phases accounting for the portions of 
the material by determining and using incident angles of ultrasound relative to the layers 
of the portions of the material. 

15 

2 1 . The computer program product of claim 2 1 wherein the instructions for 
causing the computer to determine incident angles cause the computer to determine 
incident angles for each layer of the portions of the material. 
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